（2018-09-20：Julia 1.0用にアップデート） （2018-11-01：スペクトログラム作った）

In [8]:
using Plots
pyplot();
using Images

using SampledSignals
SampledSignals.embed_javascript();

using FileIO: load, save, loadstreaming, savestreaming
import LibSndFile

using DSP
using FFTW
using Statistics

using Distributed, SharedArrays

まずは音を読み込みます。

In [9]:
snd = load("electric_bass.flac");
fs = snd.samplerate;
x = snd.data;

全体像を見てみましょう。LibSndFileの機能（というか、その土台になっているSampledSignalsの機能）を使うと、簡易的に波形表示を見ることができます。行末にセミコロンを付けなければ計算結果が表示されるのと同じようなかんじです。

In [10]:
snd
Out[10]:

SampleBuf display requires javascript

To enable for the whole notebook select "Trust Notebook" from the "File" menu. You can also trust this cell by re-running it. You may also need to re-run `using SampledSignals` if the module is not yet loaded in the Julia kernel, or `SampledSignals.embed_javascript()` if the Julia module is loaded but the javascript isn't initialized.

自分でプロットするには以下のようにします。

In [11]:
t = (0:length(x)-1) / fs;
plot(t, x,
    xlabel="Time (s)", ylabel="Amplitude",
    xticks=0:10:t[end],
    ylim=(-1, +1),
    legend=false)
Out[11]:

これは開放弦（E1、A1、D2、G2）を順番に一回ずつ演奏したもので、一つのファイルに4音が入っています。先頭にあるE1音のみを、3.85秒〜41.00秒と範囲指定して抜き出します。ただし秒数ではなくサンプル数で指定する必要があるので、秒数に標本化周波数をかけてサンプル数を計算します。また、そのままだと実数型（Float）なので、整数型（Int）に変換するためにround()を使っています。

xE1はE1音の振幅が、tE1は各サンプルの時刻が入ります。

In [12]:
xE1 = x[round(Int, 3.850*fs) : round(Int, 38.960*fs)];
tE1 = (0:length(xE1)-1) / fs;
plot(tE1, xE1,
    xlabel="Time (s)", ylabel="Amplitude",
    ylim=(-1, +1),
    legend=false)
Out[12]:

波形の詳細がどういう風になっているのか調べるために、一部分だけを拡大してみます。たとえば、x軸で2.7秒〜2.9秒の部分だけを拡大表示するときにはxlim=(2.7, 2.9)とします。

In [13]:
plot(tE1, xE1,
    xlabel="Time (s)", ylabel="Amplitude",
    xlim=(2.7, 2.9),
    ylim=(-1, +1),
    legend=false)
Out[13]:
In [14]:
yE1 = fft(xE1);
ampspecE1 = abs.(yE1[1:round(Int, length(yE1)/2+1)]);
ampspecE1[2:round(Int, length(ampspecE1)-1)] *= 2.0;

f = (0:length(ampspecE1)) / length(ampspecE1) * (fs / 2);
ampspecE1dB = 20 * log10.(ampspecE1);
plot(f, ampspecE1dB .- maximum(ampspecE1dB),
    xlim=[0, 1500], ylim=[-75, 3],
    xlabel="Frequency [Hz]", ylabel="Power [dB]",
    title="Electric Bass (E1, 41 Hz)", legend=false)
Out[14]:

倍音成分がだいたい41の倍数の周波数のところに出ていることが分かります。また、400 Hz付近と800 Hz付近に谷ができていますが、これは弾弦位置によるものです。

フレームごとに分割してスペクトル計算

In [15]:
framesize = 2^14;
framestep = 2^13;
Nmax = 2^18;

xx = zeros((length(xE1) + framesize));
xx[1:length(xE1)] = xE1;
w = hanning(framesize);

zz = SharedArray{Float64}(floor(Int, length(xx)/framestep), round(Int, max(nextpow(2, framesize), Nmax)/2+1));
@sync @distributed for n = 1:floor(Int, length(xE1)/framestep)
    xxx = xx[(n-1)*framestep+1 : (n-1)*framestep+framesize];
    z = zeros(max(nextpow(2, framesize), Nmax));
    z[1:framesize] = xxx .* w;
    y = FFTW.fft(z) / length(z);
    y[2:round(Int, length(y)/2)] *= 2.0;
    y = y[1:round(Int, length(y)/2+1)];
    zz[n, :] = abs.(y);
end
t = range(0, stop=length(xx)-1, length=size(zz, 1)) / fs;
f = range(0, stop=fs/2, length=size(zz, 2));
In [ ]:
size(zz)
In [16]:
zf_dB = 20*log10.(mean(zz, dims=1)');
plot(f, zf_dB .- maximum(zf_dB),
    xlim=(0, 1500), ylim=(-75, 3),
    xlabel="Frequency (Hz)", ylabel="Power (dB)",
    title="Electric Bass (E1, 41 Hz)", legend=false)
Out[16]:
In [17]:
zt_dB = 20*log10.(mean(zz, dims=2));
plot(t, zt_dB .- maximum(zt_dB),
    ylim=(-75, 3),
    xlabel="Time (s)", ylabel="Power (dB)",
    title="Electric Bass (E1, 41 Hz)", legend=false)
Out[17]:
In [18]:
heatmap(t, f, 20*log10.(zz' .+ eps()), xlabel="Time (s)", ylabel="Frequency (Hz)", ylim=(0,1500))
Out[18]:
In [20]:
savefig("heatmap.png")
In [ ]:
zmax = 20*log10.(maximum(zz));

anim = @animate for n=1:size(zz, 1)
    zframe = 20*log10.(zz[n, :]));
    plot(f, zframe .- zmax,
        xlim=(0, 1500), ylim=(-75, 3),
        xlabel="Frequency (Hz)", ylabel="Power (dB)",
        title="Electric Bass (E1, 41 Hz)", legend=false)
end
gif(anim, "ebass.gif", fps=15)